//本文件是普通高斯消去pthread并行化有关的各函数的合集
//注意，进行测试时可能需要分到多个文件中，由于函数来自多个工程，所以有一些变量的冲突什么的

/*
****我实现的算法有***
normal              avx             未pthread并行化
normal_8            avx_8           基础pthread并行化，无同步
normal_sem          avx_sem         pthread并行化，信号量同步
normal_semthr       avx_semthr      pthread并行化，信号量同步，三重循环
normal_barrier      avx_barrier     pthread并行化，barrier同步
normal_dy           avx_dy          pthread并行化，动态（创建）线程
以上都是循环分配，水平划分
**********************
****不同任务划分方式***
normal_sem_ver      normal_sem的水平划分变垂直划分
normal_sem_block    normal_sem的循环分配变数据块
**********************
*/

#include<iostream>
#include <stdio.h>
#include<typeinfo>
#include <stdlib.h>
#include<semaphore.h>
#include<pthread.h>
#include<tmmintrin.h>
#include<xmmintrin.h>
#include<emmintrin.h>
#include<pmmintrin.h>
#include<smmintrin.h>
#include<nmmintrin.h>
#include<immintrin.h>
#include<windows.h>
using namespace std;
int N = 100; //问题规模

#define NUM_THREADS 7
float** A = NULL;

long long head, tail, freq; //测时间

sem_t sem_main; 
sem_t sem_workstart[NUM_THREADS];
sem_t sem_workend[NUM_THREADS];

sem_t sem_leader;
sem_t sem_Division[NUM_THREADS];
sem_t sem_Elimination[NUM_THREADS];

sem_t sem_main_loop2;
sem_t sem_workerstart_loop2[NUM_THREADS];
sem_t sem_workerend_loop2[NUM_THREADS];

sem_t sem_main_block;
sem_t sem_workerstart_block[NUM_THREADS];
sem_t sem_workerend_block[NUM_THREADS];

pthread_barrier_t barrier_Division;
pthread_barrier_t barrier_Elimination;

// 参数数据结构
struct threadParam_t {
    int k;
    int t_id;
};

// 未对齐的数组的初始化
void A_init() {
    A = new float* [N];
    for (int i = 0; i < N; i++) {
        A[i] = new float[N];
    }
    for (int i = 0; i < N; i++) {
        A[i][i] = 1.0;
        for (int j = i + 1; j < N; j++) {
            A[i][j] = rand() % 1000;
        }

    }
    for (int k = 0; k < N; k++) {
        for (int i = k + 1; i < N; i++) {
            for (int j = 0; j < N; j++) {
                A[i][j] += A[k][j];
                A[i][j] = (int)A[i][j] % 1000;
            }
        }
    }
}

// 释放矩阵空间
void deleteA() {
    for (int i = 0; i < N; i++) {
        delete[] A[i];
    }
    delete A;
}

// 打印，测试
void print(float** a) {
    for (int i = 0; i < N; i++) {
        for (int j = 0; j < N; j++) {
            cout << a[i][j] << " ";
        }
        cout << endl;
    }
}

//调用并计时
void cal(void(*func)()) {
    A_init();
    QueryPerformanceCounter((LARGE_INTEGER*)&head);
    func();
    QueryPerformanceCounter((LARGE_INTEGER*)&tail);
}

//01 普通算法串行
void normal() {
    for (int k = 0; k < N; k++) {
        for (int j = k + 1; j < N; j++) {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;

        for (int i = k + 1; i < N; i++) {
            for (int j = k + 1; j < N; j++) {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;
        }
    }
}

//02 SIMD_avx优化算法
void normal_avx() {
    for (int k = 0; k < N; k++) {
        __m256 vt = _mm256_set1_ps(A[k][k]);
        int j = 0;
        for (j = k + 1; j + 8 <= N; j += 8) {
            __m256 va = _mm256_loadu_ps(&A[k][j]);
            va = _mm256_div_ps(va, vt);
            _mm256_storeu_ps(&A[k][j], va);
        }
        for (; j < N; j++) {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;
        for (int i = k + 1; i < N; i++) {
            __m256 vaik = _mm256_set1_ps(A[i][k]);
            for (j = k + 1; j + 8 <= N; j += 8) {
                __m256 vakj = _mm256_loadu_ps(&A[k][j]);
                __m256 vaij = _mm256_loadu_ps(&A[i][j]);
                __m256 vx = _mm256_mul_ps(vakj, vaik);
                vaij = _mm256_sub_ps(vaij, vx);
                _mm256_storeu_ps(&A[i][j], vaij);
            }
            for (; j < N; j++) {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0;

        }
    }
}

//03 pthread并行化，动态（创建）线程，循环水平
void* normal_dyfunc(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int k = p->k;           //消去的轮次
    int t_id = p->t_id;     //线程
    int i = k + t_id + 1;   //获取任务

    for (int j = k + 1; j < N; j++) {
        A[i][j] = A[i][j] - A[i][k] * A[k][j];
    }
    A[i][k] = 0;
    pthread_exit(NULL);
    return NULL;
}
void normal_dy() {
    for (int k = 0; k < N; k++) {
        for (int j = k + 1; j < N; j++) {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;

        int thread_cnt = N - 1 - k;
        pthread_t* handle = (pthread_t*)malloc(thread_cnt * sizeof(pthread_t));
        threadParam_t* param = (threadParam_t*)malloc(thread_cnt * sizeof(threadParam_t));

        for (int t_id = 0; t_id < thread_cnt; t_id++) {//分配任务
            param[t_id].k = k;
            param[t_id].t_id = t_id;
        }

        for (int t_id = 0; t_id < thread_cnt; t_id++) {
            pthread_create(&handle[t_id], NULL, normal_dyfunc, &param[t_id]);
        }

        for (int t_id = 0; t_id < thread_cnt; t_id++) {
            pthread_join(handle[t_id], NULL);
        }
        free(handle);
        free(param);
    }

}

//04 avx优化，pthread并行化，动态（创建）线程，循环水平
void* avx_dyfunc(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int k = p->k;           //消去的轮次
    int t_id = p->t_id;     //线程
    int i = k + t_id + 1;   //获取任务

    __m256 vaik = _mm256_set1_ps(A[i][k]);
    int j;
    for (j = k + 1; j + 8 <= N; j += 8) {
        __m256 vakj = _mm256_loadu_ps(&A[k][j]);
        __m256 vaij = _mm256_loadu_ps(&A[i][j]);
        __m256 vx = _mm256_mul_ps(vakj, vaik);
        vaij = _mm256_sub_ps(vaij, vx);
        _mm256_storeu_ps(&A[i][j], vaij);
    }
    for (; j < N; j++) {
        A[i][j] = A[i][j] - A[i][k] * A[k][j];
    }
    A[i][k] = 0;
    pthread_exit(NULL);
    return NULL;
}
void avx_dy() {
    for (int k = 0; k < N; k++) {
        __m256 vt = _mm256_set1_ps(A[k][k]);
        int j = 0;
        for (j = k + 1; j + 8 <= N; j += 8) {
            __m256 va = _mm256_loadu_ps(&A[k][j]);
            va = _mm256_div_ps(va, vt);
            _mm256_storeu_ps(&A[k][j], va);
        }
        for (; j < N; j++) {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;

        int thread_cnt = N - 1 - k;
        pthread_t* handle = (pthread_t*)malloc(thread_cnt * sizeof(pthread_t));
        threadParam_t* param = (threadParam_t*)malloc(thread_cnt * sizeof(threadParam_t));

        for (int t_id = 0; t_id < thread_cnt; t_id++) {
            param[t_id].k = k;
            param[t_id].t_id = t_id;
        }

        for (int t_id = 0; t_id < thread_cnt; t_id++) {
            pthread_create(&handle[t_id], NULL, avx_dyfunc, &param[t_id]);
        }

        for (int t_id = 0; t_id < thread_cnt; t_id++) {
            pthread_join(handle[t_id], NULL);
        }
        free(handle);
        free(param);
    }


}

//05 静态8线程，循环水平
void* normal_8func(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int k = p->k;           //消去的轮次
    int t_id = p->t_id;     //线程
    int i = k + t_id + 1;   //获取任务
    for (; i < N; i += NUM_THREADS) {
        for (int j = k + 1; j < N; j++) {
            A[i][j] = A[i][j] - A[i][k] * A[k][j];
        }
        A[i][k] = 0;
    }
    pthread_exit(NULL);
    return NULL;
}
void normal_8() {
    for (int k = 0; k < N; k++) {
        for (int j = k + 1; j < N; j++)
            A[k][j] = A[k][j] / A[k][k];
        A[k][k] = 1.0;

        int thread_cnt = NUM_THREADS;
        pthread_t* handle = (pthread_t*)malloc(thread_cnt * sizeof(pthread_t));
        threadParam_t* param = (threadParam_t*)malloc(thread_cnt * sizeof(threadParam_t));

        for (int t_id = 0; t_id < thread_cnt; t_id++) {
            param[t_id].k = k;
            param[t_id].t_id = t_id;
        }

        for (int t_id = 0; t_id < thread_cnt; t_id++)
            pthread_create(&handle[t_id], NULL, normal_8func, &param[t_id]);
        for (int t_id = 0; t_id < thread_cnt; t_id++) 
            pthread_join(handle[t_id], NULL);

        free(handle);
        free(param);
    }

}

//06 avx优化，静态8线程，循环水平
void* avx_8func(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int k = p->k;           //消去的轮次
    int t_id = p->t_id;     //线程
    int i = k + t_id + 1;   //获取任务

    for (; i < N; i += NUM_THREADS) {
        __m256 vaik = _mm256_set1_ps(A[i][k]);
        int j;
        for (j = k + 1; j + 8 <= N; j += 8) {
            __m256 vakj = _mm256_loadu_ps(&A[k][j]);
            __m256 vaij = _mm256_loadu_ps(&A[i][j]);
            __m256 vx = _mm256_mul_ps(vakj, vaik);
            vaij = _mm256_sub_ps(vaij, vx);
            _mm256_storeu_ps(&A[i][j], vaij);
        }
        for (; j < N; j++) {
            A[i][j] = A[i][j] - A[i][k] * A[k][j];
        }
        A[i][k] = 0;
    }
    pthread_exit(NULL);
    return NULL;
}
void avx_8() {
    for (int k = 0; k < N; k++) {
        __m256 vt = _mm256_set1_ps(A[k][k]);
        int j = 0;
        for (j = k + 1; j + 8 <= N; j += 8) {
            __m256 va = _mm256_loadu_ps(&A[k][j]);
            va = _mm256_div_ps(va, vt);
            _mm256_storeu_ps(&A[k][j], va);
        }
        for (; j < N; j++) {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;

        int thread_cnt = NUM_THREADS;
        pthread_t* handle = (pthread_t*)malloc(thread_cnt * sizeof(pthread_t));
        threadParam_t* param = (threadParam_t*)malloc(thread_cnt * sizeof(threadParam_t));

        for (int t_id = 0; t_id < thread_cnt; t_id++) {
            param[t_id].k = k;
            param[t_id].t_id = t_id;
        }

        for (int t_id = 0; t_id < thread_cnt; t_id++) {
            pthread_create(&handle[t_id], NULL, avx_8func, &param[t_id]);
        }

        for (int t_id = 0; t_id < thread_cnt; t_id++) {
            pthread_join(handle[t_id], NULL);
        }
        free(handle);
        free(param);
    }
}

//07 静态8线程+信号量，循环水平
void* normal_8_semfunc(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;
    for (int k = 0; k < N; k++) {
        sem_wait(&sem_workstart[t_id]);//在开始消去操作前确认主线程已完成除法操作
        //ID=t_id的线程从第k+1+t_id行开始，处理每隔7行的数据
        for (int i = k + 1 + t_id; i < N; i += NUM_THREADS) {
            for (int j = k + 1; j < N; j++)
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            A[i][k] = 0;
        }
        sem_post(&sem_main);//完成消去操作后通知主线程
        sem_wait(&sem_workend[t_id]);//完成当前消去后等待进入下一轮消去的信号量
    }
    pthread_exit(NULL);
    return NULL;
}
void normal_8_sem() {
    sem_init(&sem_main, 0, 0); //初始化信号量
    for (int i = 0; i < NUM_THREADS; i++) {
        sem_init(&sem_workend[i], 0, 0);
        sem_init(&sem_workstart[i], 0, 0);
    }
    pthread_t* handle = (pthread_t*)malloc(NUM_THREADS * sizeof(pthread_t));
    threadParam_t* param = (threadParam_t*)malloc(NUM_THREADS * sizeof(threadParam_t));
    for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
        param[t_id].t_id = t_id;
        param[t_id].k = 0;
        pthread_create(&handle[t_id], NULL, normal_8_semfunc, &param[t_id]);
    }
    for (int k = 0; k < N; k++) {
        for (int j = k + 1; j < N; j++)
            A[k][j] = A[k][j] / A[k][k];
        A[k][k] = 1.0;  //消去开始前，归一化
        for (int t_id = 0; t_id < NUM_THREADS; t_id++)//唤醒子线程开始消去
            sem_post(&sem_workstart[t_id]);
        for (int t_id = 0; t_id < NUM_THREADS; t_id++)//等待所有子线程完成消去操作
            sem_wait(&sem_main);
        for (int t_id = 0; t_id < NUM_THREADS; t_id++)//唤醒所有子线程进入下一轮消去
            sem_post(&sem_workend[t_id]);
    }
    for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
        pthread_join(handle[t_id], NULL);
    }
    sem_destroy(&sem_main);//销毁线程
    for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        sem_destroy(&sem_workstart[t_id]);
    for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        sem_destroy(&sem_workend[t_id]);
    free(handle);
    free(param);
}

//08 avx静态8线程+信号量，循环水平
void* avx_8_semfunc(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;

    for (int k = 0; k < N; k++) {
        sem_wait(&sem_workstart[t_id]);

        for (int i = k + 1 + t_id; i < N; i += NUM_THREADS) {
            __m256 vaik = _mm256_set1_ps(A[i][k]);
            int j;
            for (j = k + 1; j + 8 <= N; j += 8) {
                __m256 vakj = _mm256_loadu_ps(&A[k][j]);
                __m256 vaij = _mm256_loadu_ps(&A[i][j]);
                __m256 vx = _mm256_mul_ps(vakj, vaik);
                vaij = _mm256_sub_ps(vaij, vx);
                _mm256_storeu_ps(&A[i][j], vaij);
            }
            for (; j < N; j++) {
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            }
            A[i][k] = 0.0;
        }

        sem_post(&sem_main); //唤醒主线程
        sem_wait(&sem_workend[t_id]);  //阻塞，等待主线程唤醒进入下一轮

    }
    pthread_exit(NULL);
    return NULL;
}
void avx_8_sem() {
    sem_init(&sem_main, 0, 0);
    for (int i = 0; i < NUM_THREADS; i++) {
        sem_init(&sem_workend[i], 0, 0);
        sem_init(&sem_workstart[i], 0, 0);
    }
    pthread_t* handle = (pthread_t*)malloc(NUM_THREADS * sizeof(pthread_t));
    threadParam_t* param = (threadParam_t*)malloc(NUM_THREADS * sizeof(threadParam_t));
    for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
        param[t_id].t_id = t_id;
        param[t_id].k = 0;
        pthread_create(&handle[t_id], NULL, avx_8_semfunc, &param[t_id]);

    }

    for (int k = 0; k < N; k++) {

        __m256 vt = _mm256_set1_ps(A[k][k]);
        int j = 0;
        for (j = k + 1; j + 8 <= N; j += 8) {
            __m256 va = _mm256_loadu_ps(&A[k][j]);
            va = _mm256_div_ps(va, vt);
            _mm256_storeu_ps(&A[k][j], va);
        }
        for (; j < N; j++) {
            A[k][j] = A[k][j] / A[k][k];
        }
        A[k][k] = 1.0;

        for (int t_id = 0; t_id < NUM_THREADS; t_id++) {  //唤起子线程
            sem_post(&sem_workstart[t_id]);
        }

        for (int t_id = 0; t_id < NUM_THREADS; t_id++) {  //主线程睡眠
            sem_wait(&sem_main);
        }

        for (int t_id = 0; t_id < NUM_THREADS; t_id++) {  //再次唤起工作线程，进入下一轮消去
            sem_post(&sem_workend[t_id]);
        }

    }
    for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
        pthread_join(handle[t_id], NULL);
    }
    sem_destroy(&sem_main);    //销毁线程
    for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        sem_destroy(&sem_workstart[t_id]);
    for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        sem_destroy(&sem_workend[t_id]);

    free(handle);
    free(param);

}

//09 静态8线程+信号量三重循环，循环水平
void* normal_8_semthrfunc(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;
    for (int k = 0; k < N; k++) { // 逐轮次进行消去操作
        if (t_id == 0) { // 线程0进行除法操作
            for (int j = k + 1; j < N; j++)
                A[k][j] = A[k][j] / A[k][k];
            A[k][k] = 1.0;
        }
        else { // 其余线程等待线程0完成除法操作
            sem_wait(&sem_Division[t_id - 1]);
        }
        if (t_id == 0) { // 线程0唤醒其他线程进行消去操作
            for (int i = 0; i < NUM_THREADS - 1; i++)
                sem_post(&sem_Division[i]);
        }
        // 各线程进行消去操作
        for (int i = k + 1 + t_id; i < N; i += NUM_THREADS) {
            for (int j = k + 1; j < N; j++)
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            A[i][k] = 0;
        }
        if (t_id == 0) { // 线程0等待其他线程完成消去操作
            for (int i = 0; i < NUM_THREADS - 1; i++)
                sem_wait(&sem_leader);
            for (int i = 0; i < NUM_THREADS - 1; i++)
                sem_post(&sem_Elimination[i]);
        }
        else { // 其他线程通知线程0消去操作完成
            sem_post(&sem_leader);
            sem_wait(&sem_Elimination[t_id - 1]);
        }
    }
    pthread_exit(NULL);
    return NULL;
}
void normal_8_semthr() {
    sem_init(&sem_leader, 0, 0); // 初始化领导信号量
    for (int i = 0; i < NUM_THREADS; i++) {
        sem_init(&sem_Division[i], 0, 0); // 初始化除法信号量
        sem_init(&sem_Elimination[i], 0, 0); // 初始化消去信号量
    }
    pthread_t* handle = (pthread_t*)malloc(NUM_THREADS * sizeof(pthread_t));
    threadParam_t* param = (threadParam_t*)malloc(NUM_THREADS * sizeof(threadParam_t));
    for (int t_id = 0; t_id < NUM_THREADS; t_id++) { // 创建线程
        param[t_id].t_id = t_id;
        param[t_id].k = 0;
        pthread_create(&handle[t_id], NULL, normal_8_semthrfunc, &param[t_id]);
    }
    for (int t_id = 0; t_id < NUM_THREADS; t_id++) // 等待所有线程结束
        pthread_join(handle[t_id], NULL);
    sem_destroy(&sem_leader); // 销毁领导信号量
    for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
        sem_destroy(&sem_Division[t_id]); // 销毁除法信号量
        sem_destroy(&sem_Elimination[t_id]); // 销毁消去信号量
    }
    free(handle); // 释放内存
    free(param);
}

//10 avx静态8线程+信号量三重循环，循环水平
void* avx_8_semthrfunc(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;

    for (int k = 0; k < N; k++) {
        if (t_id == 0) {
            __m256 vt = _mm256_set1_ps(A[k][k]);
            int j = 0;
            for (j = k + 1; j + 8 <= N; j += 8) {
                __m256 va = _mm256_loadu_ps(&A[k][j]);
                va = _mm256_div_ps(va, vt);
                _mm256_storeu_ps(&A[k][j], va);
            }
            for (; j < N; j++)
                A[k][j] = A[k][j] / A[k][k];
            A[k][k] = 1.0;
        }
        else
            sem_wait(&sem_Division[t_id - 1]);
        if (t_id == 0) {
            for (int i = 0; i < NUM_THREADS - 1; i++)
                sem_post(&sem_Division[i]);
        }
        for (int i = k + 1 + t_id; i < N; i += NUM_THREADS) {
            __m256 vaik = _mm256_set1_ps(A[i][k]);
            int j;
            for (j = k + 1; j + 8 <= N; j += 8) {
                __m256 vakj = _mm256_loadu_ps(&A[k][j]);
                __m256 vaij = _mm256_loadu_ps(&A[i][j]);
                __m256 vx = _mm256_mul_ps(vakj, vaik);
                vaij = _mm256_sub_ps(vaij, vx);
                _mm256_storeu_ps(&A[i][j], vaij);
            }
            for (; j < N; j++)
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            A[i][k] = 0.0;
        }
        if (t_id == 0) {
            for (int i = 0; i < NUM_THREADS - 1; i++)
                sem_wait(&sem_leader);
            for (int i = 0; i < NUM_THREADS - 1; i++)
                sem_post(&sem_Elimination[i]);
        }
        else {
            sem_post(&sem_leader);
            sem_wait(&sem_Elimination[t_id - 1]);
        }
    }

    pthread_exit(NULL);
    return NULL;
}
void avx_8_semthr() {
    sem_init(&sem_leader, 0, 0);
    for (int i = 0; i < NUM_THREADS; i++) {
        sem_init(&sem_Division[i], 0, 0);
        sem_init(&sem_Elimination[i], 0, 0);
    }
    pthread_t* handle = (pthread_t*)malloc(NUM_THREADS * sizeof(pthread_t));
    threadParam_t* param = (threadParam_t*)malloc(NUM_THREADS * sizeof(threadParam_t));

    for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
        param[t_id].t_id = t_id;
        param[t_id].k = 0;
        pthread_create(&handle[t_id], NULL, avx_8_semthrfunc, &param[t_id]);

    }
    for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
        pthread_join(handle[t_id], NULL);
    }
    sem_destroy(&sem_main);
    for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        sem_destroy(&sem_workstart[t_id]);
    for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        sem_destroy(&sem_workend[t_id]);

    free(handle);
    free(param);
}

//11 静态barrier，循环水平
void* normal_barrierfunc(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;
    for (int k = 0; k < N; k++) {
        if (t_id == 0) {  // 线程0负责除法操作
            for (int j = k + 1; j < N; j++)
                A[k][j] = A[k][j] / A[k][k];
            A[k][k] = 1.0;
        }
        // 第一个同步点，等待所有线程完成除法操作
        pthread_barrier_wait(&barrier_Division);
        // 消去操作并行进行
        for (int i = k + 1 + t_id; i < N; i += NUM_THREADS) {
            for (int j = k + 1; j < N; j++)
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            A[i][k] = 0;
        }
        // 第二个同步点，等待所有线程完成消去操作
        pthread_barrier_wait(&barrier_Elimination);
    }
    pthread_exit(NULL);
    return NULL;
}
void normal_barrier() {
    // 初始化屏障变量
    pthread_barrier_init(&barrier_Division, NULL, NUM_THREADS);
    pthread_barrier_init(&barrier_Elimination, NULL, NUM_THREADS);
    pthread_t* handle = (pthread_t*)malloc(NUM_THREADS * sizeof(pthread_t));
    threadParam_t* param = (threadParam_t*)malloc(NUM_THREADS * sizeof(threadParam_t));
    // 创建线程
    for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
        param[t_id].t_id = t_id;
        param[t_id].k = 0;
        pthread_create(&handle[t_id], NULL, normal_barrierfunc, &param[t_id]);
    }
    // 等待所有线程完成
    for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        pthread_join(handle[t_id], NULL);
    // 销毁屏障变量
    pthread_barrier_destroy(&barrier_Division);
    pthread_barrier_destroy(&barrier_Elimination);
    free(handle);
    free(param);
}

//12 avx静态barrier，循环水平
void* avx_barrierfunc(void* param) {
    threadParam_t* p = (threadParam_t*)param;
    int t_id = p->t_id;
    for (int k = 0; k < N; k++) {
        if (t_id == 0) {
            __m256 vt = _mm256_set1_ps(A[k][k]);
            int j = 0;
            for (j = k + 1; j + 8 <= N; j += 8) {
                __m256 va = _mm256_loadu_ps(&A[k][j]);
                va = _mm256_div_ps(va, vt);
                _mm256_storeu_ps(&A[k][j], va);
            }
            for (; j < N; j++)
                A[k][j] = A[k][j] / A[k][k];
            A[k][k] = 1.0;
        }
        pthread_barrier_wait(&barrier_Division);
        for (int i = k + 1 + t_id; i < N; i += NUM_THREADS) {
            __m256 vaik = _mm256_set1_ps(A[i][k]);
            int j;
            for (j = k + 1; j + 8 <= N; j += 8) {
                __m256 vakj = _mm256_loadu_ps(&A[k][j]);
                __m256 vaij = _mm256_loadu_ps(&A[i][j]);
                __m256 vx = _mm256_mul_ps(vakj, vaik);
                vaij = _mm256_sub_ps(vaij, vx);
                _mm256_storeu_ps(&A[i][j], vaij);
            }
            for (; j < N; j++)
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            A[i][k] = 0.0;
        }
        pthread_barrier_wait(&barrier_Elimination);
    }
    pthread_exit(NULL);
    return NULL;
}
void avx_barrier() {
    pthread_barrier_init(&barrier_Division, NULL, NUM_THREADS);
    pthread_barrier_init(&barrier_Elimination, NULL, NUM_THREADS);
    pthread_t* handle = (pthread_t*)malloc(NUM_THREADS * sizeof(pthread_t));
    threadParam_t* param = (threadParam_t*)malloc(NUM_THREADS * sizeof(threadParam_t));
    for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
        param[t_id].t_id = t_id;
        param[t_id].k = 0;
        pthread_create(&handle[t_id], NULL, avx_barrierfunc, &param[t_id]);
    }
    for (int t_id = 0; t_id < NUM_THREADS; t_id++)
        pthread_join(handle[t_id], NULL);
    pthread_barrier_destroy(&barrier_Division);
    pthread_barrier_destroy(&barrier_Elimination);
    free(handle);
    free(param);
}

//13 normal_sem_ver
void* normal_sem_verfunc(void* param) {
	threadParam_t2* p = (threadParam_t2*)param;
	int t_id = p->t_id;
	for (int k = 0; k < N; ++k) {
		sem_wait(&sem_workerstart_loop2[t_id]);//阻塞，等待主线完成除法操作（操作自己专属的信号量）
		//循环任务划分	垂直划分
		for (int i = k + 1; i < N; i ++) {
			for (int j = k + 1 + t_id; j < N; j+=NUM_THREADS)
				A[i][j] = A[i][j] - A[i][k] * A[k][j];
			A[i][k] = 0.0;
		}
		sem_post(&sem_main_loop2);//唤醒主线程
		sem_wait(&sem_workerend_loop2[t_id]);//阻塞，等待主线程唤醒进入下一轮
	}
	pthread_exit(NULL);
	return NULL;
}
void normal_sem_ver() {
	//初始化信号量
	sem_init(&sem_main_loop2, 0, 0);
	for (int i = 0; i < NUM_THREADS; ++i) {
		sem_init(&sem_workerstart_loop2[i], 0, 0);
		sem_init(&sem_workerend_loop2[i], 0, 0);
	}
	//创建线程
	pthread_t handles[NUM_THREADS];
	threadParam_t2 param[NUM_THREADS];
	for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
		param[t_id].t_id = t_id;
		pthread_create(&handles[t_id], NULL, normal_sem_verfunc, (void*)&param[t_id]);
	}
	for (int k = 0; k < N; k++) {
		//主线程做除法操作
		for (int j = k + 1; j < N; j++)
			A[k][j] = A[k][j] / A[k][k];
		A[k][k] = 1.0;
		//开始唤醒工作线程
		for (int t_id = 0; t_id < NUM_THREADS; t_id++)
			sem_post(&sem_workerstart_loop2[t_id]);
		//主线程睡眠（等待所有的工作线程完成此轮消去任务）
		for (int t_id = 0; t_id < NUM_THREADS; t_id++)
			sem_wait(&sem_main_loop2);
		//主线程再次唤醒工作线程进入下一轮的消去任务
		for (int t_id = 0; t_id < NUM_THREADS; t_id++)
			sem_post(&sem_workerend_loop2[t_id]);
	}
	for (int t_id = 0; t_id < NUM_THREADS; t_id++)
		pthread_join(handles[t_id], NULL);
	//销毁所有信号量
	sem_destroy(&sem_main_loop2);
	sem_destroy(&sem_workerstart_loop2[NUM_THREADS]);
	sem_destroy(&sem_workerend_loop2[NUM_THREADS]);
    free(handles);
    free(param);
}

//14 normal_sem_block
void* normal_sem_blockfunc(void* param) {
	threadParam_t2* p = (threadParam_t2*)param;
	int t_id = p->t_id;
	for (int k = 0; k < N; ++k) {
		sem_wait(&sem_workerstart_block[t_id]);//阻塞，等待主线完成除法操作（操作自己专属的信号量）
		//任务按块划分
		int q = (N - k - 1) / (NUM_THREADS+1);
		int w = k + t_id * q + 1; //获取自己的计算任务
		for (int i = w; i < w + q; i++)
		{
			for (int j = k + 1; j < N; ++j)
			{
				A[i][j] = A[i][j] - A[i][k] * A[k][j];
			}
			A[i][k] = 0;
		}
		sem_post(&sem_main_block);//唤醒主线程
		sem_wait(&sem_workerend_block[t_id]);//阻塞，等待主线程唤醒进入下一轮
	}
	pthread_exit(NULL);
	return NULL;
}
void normal_sem_block() {
	//初始化信号量
	sem_init(&sem_main_block, 0, 0);
	for(int i = 0; i < NUM_THREADS; ++i) {
		sem_init(&sem_workerstart_block[i], 0, 0);
		sem_init(&sem_workerend_block[i], 0, 0);
	}
	//创建线程
	pthread_t handles[NUM_THREADS];
	threadParam_t2 param[NUM_THREADS];
	for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
		param[t_id].t_id = t_id;
		pthread_create(&handles[t_id], NULL, normal_sem_blockfunc, (void*)&param[t_id]);
	}
	for (int k = 0; k < N; k++) {
		//主线程做除法操作
		for (int j = k + 1; j < N; j++)
			A[k][j] = A[k][j] / A[k][k];
		A[k][k] = 1.0;
		//开始唤醒工作线程
		for (int t_id = 0; t_id < NUM_THREADS; t_id++)
			sem_post(&sem_workerstart_block[t_id]);
        for (int i = k + (N - k - 1) / (NUM_THREADS+1) * NUM_THREADS + 1; i < N; i++) {
            for (int j = k + 1; j < N; ++j)
                A[i][j] = A[i][j] - A[i][k] * A[k][j];
            A[i][k] = 0.0;
        }
		//主线程睡眠（等待所有的工作线程完成此轮消去任务）
		for (int t_id = 0; t_id < NUM_THREADS; t_id++)
			sem_wait(&sem_main_block);
		//主线程再次唤醒工作线程进入下一轮的消去任务
		for (int t_id = 0; t_id < NUM_THREADS; t_id++)
			sem_post(&sem_workerend_block[t_id]);
	}
	for (int t_id = 0; t_id < NUM_THREADS; t_id++) {
		pthread_join(handles[t_id], NULL);
	}
	//销毁所有信号量
	sem_destroy(&sem_main_block);
	sem_destroy(&sem_workerstart_block[NUM_THREADS]);
	sem_destroy(&sem_workerend_block[NUM_THREADS]);
    free(handles);
    free(param);
}


int main()
{
    QueryPerformanceFrequency((LARGE_INTEGER*)&freq);
    cal(normal);
    cout << "normal:" << (tail - head) * 1000 / freq << "ms" << endl;
    deleteA();
    return 0;
}
